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Monte Carlo simulations have been carried out for a system of monomers on square lattices that, by decreasing 
temperature or increasing density, polymerize reversibly into chains with two allowed directions and, at the same 
time, undergo a continuous isotropic-nematic (IN) transition. The results show that the self-assembly process 
affects the nature of the transition. Thus, the calculation of the critical exponents and the behavior of Binder 
cumulants indicate that the universality class of the IN transition changes from two-dimensional Ising-type for 
monodisperse rods without self-assembly to ^ = 1 Potts-type for self-assembled rods. 
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Self-assembly is a challenging field of research, driven 
principally by the desire to design new materials. Moreover, 
self-assembly is used permanently in biological systems to 
construct supramolecular structures such as virus capsids, fil- 
aments, and many others large molecular complexes. So, un- 
derstanding the rules of self-assembly has important applica- 
tions to both materials science and biology [1]. 

On the other hand, the isotropic-nematic (IN) transition in 
solutions of rodlike particles has been attracting a great deal 
of interest since long ago. A seminal contribution to this sub- 
ject was made by Onsager [2] with his paper on the IN tran- 
sitions of infinitely thin rods. This theory shows that parti- 
cles interacting with only excluded volume interaction may 
exhibit a rich phase diagram, despite the absence of any at- 
traction. Later, computer simulations of hard ellipses of finite 
length [3] confirmed the Onsager prediction that particle sha- 
pe anisotropy can be a sufficient condition to induce the long- 
range orientational order found in nematic liquid crystals. 

In contrast to ordinary liquid crystal, many rodlike bio- 
logical polymers are formed by monomers reversibly self- 
assembling into chains of arbitrary length so that these sys- 
tems exhibit a broad equilibrium distribution of filament 
lengths. An experimental contribution to the study of these 
systems has been presented by Viamontes et al. [4J. The au- 
thors reported a continuous IN transition for solutions of long 
F-actin (average filament length longer than 2/im) and showed 
the existence of a first-order phase separation for solutions of 
F-actin with average filament length shorter than 2jj.m. These 
findings contradict what is generally accepted in the litera- 
ture: in three dimensions, the IN transition is typically first 
order. On the other hand, in two dimensions both continu- 
ous Istl and first-order IN transitions can occur. Here, we 
consider a self-assembled two-dimensional (2D) system that 
undergo a IN transition, which is expected to be a continuous 
phase transition ff]. 

As mentioned above, the self-assembled system is intrinsi- 
cally polydisperse. While being able to solve explicitly only 
the monodisperse case, Onsager |2] already outlined the pos- 
sible extension of the theory to polydisperse systems. In this 



line of work, a detailed investigation of the effects of full 
length polydispersity, i.e., of a continuous distribution of rod 
lengths, on the Onsager theory has been recently developed by 
Speranza and Sollich ["s*]. Another approach to the problem of 
monodisperse rodlike mixtures has been proposed by Zwanzig 
|9]. The Zwanzig model has been also extended to polydis- 
perse systems flO], providing thus a useful starting point for 
understanding the effects of polydispersity on the phase be- 
havior of hard rod systems. However, a complete description 
of a system of self-assembled rods should consider not only 
the effects of polydispersity, but also the influence of the poly- 
merization process. 

In this context, we focus on a system composed of 
monomers with two attractive (sticky) poles that polymerize 
reversibly into polydisperse chains and, at the same time, un- 
dergo a continuous phase transition. So, the interplay between 
the self-assembly process and the nematic ordering is a dis- 
tinctive characteristic of these systems. 

The same system has been recently considered by Tavares 
et al. [7]. Using an approach in the spirit of the Zwanzig 
model, the authors studied the IN transition occurring in a 
two-dimensional system of self-assembled rigid rods. The ob- 
tained results revealed that nematic ordering enhances bond- 
ing. In addition, the average rod length was described quan- 
titatively in both phases, while the location of the ordering 
transition, which was found to be continuous, was predicted 
semiquantitatively by the theory. 

Despite these interesting results there is an open question to 
be answered: "what type of phase transition is it?" Tavares et 
al. ^ assumed as working hypothesis that the nature of the 
IN transition remains unchanged with respect to the case of 
monodisperse rigid rods on square lattices, where the transi- 
tion is in the 2D Ising universality class lfTll[l2ll . In this con- 
text, the confirmation (or not) of this hypothesis is not only 
important to resolve the universality class of the IN transi- 
tion occuiTing in a system of self-assembled rods, but also 
to shed light on our understanding of the effect of the self- 
assembly process on the nature of the transition. The objective 
of this Rapid Communication is to provide a thorough study in 
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this direction. For this purpose, extensive Monte Carlo (MC) 
simulations supplemented by analysis using finite-size scaling 
(FSS) theory Ml 311 have been carried out to study the critical 
behavior in a system of self-assembled rigid rods deposited on 
square lattices with two allowed directions. The calculations 
were developed at constant temperature and different densi- 
ties, thus allowing a direct comparison with previous results 
for l ong monodisperse rigid rods on two-dimensional lattices 
111 l[ll2ll . Then, the conventional normalized scaling variable 
e = T /Tc — 1 was replaced by e = — 1, where T , Tc, 
9, and Oc represent temperature, critical temperature, density, 
and critical density, respectively. A nematic phase, character- 
ized by a big domain of parallel self-assembled rigid rods, is 
separated from the disordered state by a continuous IN transi- 
tion occurring at a finite critical density. The results show that 
the self-assembly process affects the nature of the transition. 
Thus, the determination of the critical exponents indicates that 
the universality class of the IN transition changes from 2D 
Ising-type for monodisperse rods without self-assembly lll2ll 
to <7 = 1 Potts-type for self-assembled rods. 

In a recent paper, Fischer and Vink lll4ll indicated that the 
transition studied in Refs. Ill llll2ll corresponds to a liquid-gas 



transition rather than IN. This interpretation is consistent with 
the 2D-Ising critical behavior observed for monodisperse rigid 
rods on square lattices. However, as mentioned in the previ- 
ous paragraph, the universality class of self-assembled rods 
changes from that of the 2D Ising model and the result in Ref. 
lfr4ll is not generalizable to the system studied here. Accord- 
ingly, we will continue using the term "IN phase transition" 
as in the previous work of Tavares et al. ffT] . 

As in Ref. 101, we consider a system of self-assembled rods 
with a discrete number of orientations in 2D. We assume that 
the substrate is represented by a square lattice of M ^ Lx L 
sites with periodic boundary conditions. particles are ad- 
sorbed on the substrate with two possible orientations along 
the principal axis of the lattice. These particles interact with 
nearest neighbors (NN) through anisotropic attractive interac- 
tions. Then the adsorbed phase is characterized by the Hamil- 
tonian H = Y,{i,j)^iiCiC j, where indicates a sum over 
NN sites; Wjj represents the NN lateral interaction, which is 
Wij = —w if two neighboring particles / and j are aligned with 
each other and with the intermolecular vector, and is Wij — 
otherwise; and c, is the occupation variable with c, = if the 
site / is empty, and c,- = 1 if the site / is occupied. 

A cluster or uninterrupted sequence of bonded particles is 
a self-assembled rod. At fixed temperature, the average rod 
length increases as the density increases and thepolydisperse 
rods will undergo a nematic ordering transition 17|] . In order to 
follow the formation of the nematic phase from the isotropic 
phase, we use the order parameter defined in Ref. ITj], which 
can be written as 5 = \Nij -Nv\ / {Nh +Nv), where Nh{Nv) is 
the number of particles in clusters aligned along the horizontal 
(vertical) direction. When the system is disordered, all orien- 
tations are equivalents and 8 is zero. In the critical regime, 
the particles align along one direction and 5 is different from 
zero. In other words, the IN phase transition is accomplished 
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FIG. 1: Size dependence of the order parameter (5) as a function 
of density (0). Inset: size dependence of the susceptibility {%) as a 
function of density {&). 
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FIG. 2: Curves of Ul{0) vs 6 for square lattices of different sizes. 
From their intersections one obtained 6c. In the inset, the data are 
plotted over a wider range of densities. 



by a breaking of the orientational symmetry and 5 appears as 
a proper order parameter to elucidate this phenomenon. 

The problem has been studied by canonical Monte Carlo 
simulations using an vacancy -particle-exchange Kawasaki dy- 
namics ifisll and Metropolis acceptance probability lfl6ll . Typ- 
ically, the equilibrium state can be well reproduced after dis- 
carding the first 5 x 10^ Monte Carlo steps (MCS). Then, the 
next 6 X 10^ MCS are used to compute averages. 

In our Monte Carlo simulations, we set the temperature T, 
varied the density 9 — N/M, and monitored the order param- 
eter 8, which can be calculated as simple average. The quan- 
tities related with the order parameter, such as the susceptibil- 
ity and the reduced fourth-order cumulant Ul introduced 
by Binder (l^ were calculated as ;^ = I^ii^^) ~ i^)^] 
Ul = I — (5^) / [3(5^)^] , where (• • •) means the average over 
the MC simulation runs. In addition, in order to discuss the 
nature of the phase transition, the fourth-order energy cumu- 
lant Ue was obtained as f/£ = 1 - (//"*) / [3{H^-f] . 

The critical behavior of the present model has been inves- 
tigated by means of the computational scheme described in 
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FIG. 3: (a) Log-log plot of the size dependence of the maximum val- 
ues of derivatives of various thermodynamic quantities used to de- 
termine V. (b) Log-log plot of the size dependence of the maximum 
value of the susceptibility, the point of inflection of the order param- 
eter and the maximum value of the derivative of the order parameter 
used to determine yand j8, respectively. 
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FIG. 4: Data collapsing of the order parameter, SL^I^ vs |e|L'/*'. 
The plot was made using 8c = 0.524 and the exact percolation expo- 
nents V = 4/3 and j3 = 5/36. 
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FIG. 5: Data collapsing of the cumulant, Ui vs eL'/*', and of the 
susceptibility, XL^'^I^ vs eL'/*' (inset). The plots were made using 
Qc = 0.524 and exact percolation exponents v = 4/3 and J3 = 43/18. 



the previous paragraphs and finite-size scaHng analysis. The 
FSS theory implies the following behavior of 5, %, and Ul 
at criticality: 5 = L-'5/*'5(L'/*'e), % = L''/*'^(L'/*'e), and 
Vl = ULiL^l^e) for L ^ oo, e such that L'/*'£= finite, 
where e = / 0^ — 1 . Here /3 , y, and v are the standard critical 
exponents of the order parameter, susceptibility, and correla- 
tion length, respectively. 8,% ™d Ul are scaling functions for 
the respective quantities. 

The phase diagram of the system under study has been re- 
cently reported by Tavares et al. The authors showed 
that the critical density, at which the IN transition occurs, in- 
creases monotonically as /w is increased. Thus, the ne- 
matic phase is stable at low temperatures and high densities 
[see Fig. 1(a) in Ref. In addition, Tavares et al. found 
strong numerical evidence that the IN transition is a continu- 
ous phase transition. However, the authors were not able to 
determine the critical quantities characterizing the universal- 
ity class of the mentioned transition. In the following, we try 
to resolve this problem. 

In our study and based on the phase diagram given in Ref. 
10], we set the lateral interaction to w = AksT. With this 
value of w, it is expected the appearance of a nematic phase 
at intermediate densities. Accordingly, the density was var- 
ied between 0.4 and 0.6. For each value of 0, the effect of 
finite size was investigated by examining square lattices with 
L = 60, 80, 100, and 120. 

We start with the calculation of the order parameter (Fig. 
1), susceptibility (inset in Fig. 1), and cumulant (Fig. 2) plot- 
ted versus Q for several lattice sizes. In the vicinity of the 
critical point, cumulants show a strong dependence on the sys- 
tem size. However, at the critical point the cumulants adopt 
a nontrivial value U*; irrespective of system sizes in the scal- 
ing limit. Thus, plotting Ul{9) for different linear dimensions 
yields an intersection point U*, which gives an accurate es- 
timation of the critical density in the infinite system and al- 
lows us to make a preliminary identification of the universal- 
ity class of the transition [13]. In this case, the values ob- 
tained for the critical density and the intersection point of the 
cumulants were 9c = 0.524(4) and U* = 0.639(3), respec- 
tively. This fixed value of the cumulants has changed from 
that obtained for monodisperse rigid rods on square lattices 
[U* = 0.615(5)], which may be taken as a first indication that 
the universality class of the present model is different from 
the well-known 2D Ising-type for monodisperse rods [12]. In 
the lower-right inset, the data are plotted over a wider range 
of densities. As can be seen, the curves exhibit the typical 
behavior of Ul in the presence of a continuous phase transi- 
tion. Namely, the order-parameter cumulant shows a smooth 
increase from to 2/3 instead of the characteristic deep (neg- 
ative) minimum, as in a first-order phase transition [ 1 3i] . 

In order to discard the possibility that the phase transition 
is a first-order one, the energy cumulants have been measured 
for different lattice sizes ranging between L — 2Q and L = 80. 
As is well known, the finite-size analysis of Ue is a simple 
and direct way to determine the order of a phase transition 
1 13|] . Our results for Ue show a dip close to the critical density 
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Oc = 0.524 for all system sizes, but this minimum scales to 2/3 
in the thermodynamic limit as can be seen in the upper-left 
inset of Fig. 2. These results exclude a first-order transition, 
confirming the predictions by Tavares et al. |7]. 

Next, the critical exponents will be calculated. As stated in 
Ref. IitIi . the critical exponent v can be obtained by consider- 
ing the scaling behavior of certain thermodynamic derivatives 
with respect to the density 9, for example, the derivative of 
the cumulant and the logarithmic derivatives of (5) and (8^). 
In Fig. 3(a), we plot the maximum value of these derivatives 
as a function of system size on a log-log scale. The results for 
1/v from these fits are given in Fig. 3(a). Combining these 
three estimates, we obtain v = 1.33(1). Once we know v, the 
exponent 7 can be determined by scaling the maximum value 
of the susceptibility [17]. Our data for max are shown in Fig. 
3(b). The value obtained for 7 is 7=2.36(4). 

On the other hand, the standard way to extract the exponent 
ratio j3 / V is to study the scaling behavior of (5) at the point of 
inflection ((5)|;i^f), i.e., at the point where d{d)/d9 is maxi- 
mal. The scaling of (5) is shown in Fig. 3(b). The linear fit 
through all data points gives jS'^^^linf) =0.139(12). In the case 
of d{d)/dO\ 

max [see Fig. 3(b)], the value obtained from the 
fit is l}('HS)/ii0U^) ^ 0.138(3). Combining the two estimates, 
we obtain the final value (5 = 0.139(6). 

The values calculated for v, )3 and 7 clearly indicate that 
the IN phase transition belongs to the q = I Potts universality 
class (ordinary percolation) 11811 . This finding is also consis- 
tent with the crossing point of the cumulants shown in Fig. 2, 
which is in agreement with recent simulations on 2D site per- 
colation by Vink lfl9ll . where a value of U* « 0.638 has been 
obtained. 

The universality observed in this Rapid Communication 
can be interpreted by analyzing the connection between the 
thermal phase transition (IN phase transition) occurring in the 
system and the behavior of the clusters of aligned monomers. 
In fact, preliminary results (based on the observation of the 
adsorbed state at critical regime) suggest that, at intermedi- 
ate concentrations, the appearance of nematic order at criti- 
cal density is accompanied by the simultaneous presence of a 
cluster of aligned monomers connecting the extremes of the 
lattice. In these conditions, an exact correspondence is ex- 
pected between ( 1 ) the nematic order parameter and the infi- 
nite cluster; (2) the susceptibility and the mean cluster size; 
and (3) the IN phase transition and the percolation phase tran- 
sition. A similar behavior has been observed for the site- 
diluted quenched Ising model at the percolation threshold 

The results are also consistent with previous research on 
polymer systems in 2D such as polymer networks 12111 and 
branched polymers li2Zl . which were shown to be in the uni- 
versality class of ordinary percolation. Furthermore, the study 
of linear segments of size k and A:-mers of different structures 
and forms deposited on 2D regular lattices has demonstrated 
that the system, in all cases, belongs to the random percolation 
universality class lEsll . 

The scaling behavior can be further tested by plotting 



(5)l'5/^ vs |e|Ll/^ ;tL-l'/^ vs eL^/\ and U vs |e|Li/^ and 
looking for data collapsing. Using 9^ = 0.524, and the ex- 
act values of the critical exponents of the ordinary percolation 
V = 4/3, j3 = 5/36, and 7 = 43/18, we obtain an excellent 
scaling collapse as it is shown in Figs. 4 and 5. This study 
leads to independent controls and consistency checks of the 
values of all the critical exponents. 

In summary, we have used Monte Carlo simulations and 
finite-size scaling theory to resolve the nature and universality 
class of the IN phase transition occurring in a model of self- 
assembled rigid rods. The existence of a continuous phase 
transition was confirmed. In addition, as was evident from 
our results, the self-assembly process affects the universality 
of the IN transition. Thus, the accurate determination of the 
critical exponents along with the behavior of Binder cumu- 
lants revealed that the universality class of the IN transition 
changes from 2D Ising-type for monodisperse rods without 
self-assembly to <7 = 1 Potts-type for self-assembled rods. 
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